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ABSTRACT 

The sunspot drawings of Johann Staudacher of 1749-1799 were used to determine the solar differential rotation in that period. These 
drawings of the full disk lack any indication of their orientation. We used a Bayesian estimator to obtain the position angles of the 
drawings, the corresponding heliographic spot positions, a time offset between the drawings and the differential rotation parameter 
SQ, assuming the equatorial rotation period is the same as today. The drawings are grouped in pairs, and the resulting marginal 
distributions for SQ were multiplied. We obtain SQ = -0.048 ± 0.025 d _1 (-2.75 7d) for the entire period. There is no significant 
difference to the value of the present Sun. We find an (insignificant) indication for a change of SQ throughout the observing period 
from strong differential rotation, SQ ~ —0.07 d , to weaker differential rotation, SQ ~ —0.04 d _1 . 
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1. Introduction 

The solar activity cycle is the result of an oscillatory magnetic 
dynamo in the interior of the Sun. Practically all available dy- 
namo models for the Sun require a differential rotation for part 
of the amplification of the magnetic fields. Today's solar differ- 
ential rotation is known from helioseismology with fairly high 
precision down to the bottom of the convection zone. 

The rotation period of the pole is considerably longer than 
the one at the equator. Balthasar et al. (1986) determined the 
differential rotation from the rotational motion of sunspots. 
Although sunspots provide no information on the rotational pe- 
riod at the pole, we can extrapolate certain profiles of the rota- 
tion to the pole and obtain a polar rotation period, which is 25% 
longer than the equatorial one. While the profile in Balthasar 
et al. (1986) runs from about 24.7 d at the equator to 31 d at 
the pole, the rotation profile obtained from helioseismology runs 
from about 25.2 d to about 38 d (Korzennik & Eff-Darwich 
2011). It is a typical feature of solar activity that sunspot rota- 
tion periods are slightly shorter than the rotation period of the 
photosphere (Benevolenskaya et al. 1999). 

Several simulations of the dynamo process included the 
modulation of the differential rotation as a source of the solar 
cycle variability. Many global dynamo solutions are kinematic 
solution, which do not integrate an equation of motion over time. 
Small-scale motions are comprised in a source term that usu- 
ally contains a nonlinearity that mimicks the back-reaction of 
the magnetic fields on the flows (a-effect and a-quenching). 
Typically, the source term is significantly diminished for mag- 
netic fields B stronger than the equipartition field strength de- 
fined by B cq = ytfiopu', where po and p are the magnetic 
permeability and the density, respectively, and u' is the typical 
turbulent velocity. It is, however, not only the feedback to the 
small-scale motions but also the back-reaction on the large-scale 
flows, in particular the differential rotation, which introduces a 
non-linearity to the dynamo system. 

The fields generated in the dynamo process exert Lorentz 
forces, which modulate the differential rotation. This was stud- 
ied by Malkus & Proctor (1975) who found highly nonlinear 



behaviour and complex time-series. These authors applied the 
Lorentz force from the large-scale magnetic field to the large- 
scale motion (differential rotation). In a refined attempt, Kiiker et 
al. (1999) also included the back-reaction of the generated large- 
scale field on the turbulent redistribution of angular momentum 
(A-effect and A-quenching). They found time-variations of the 
dynamo solution that are very reminiscent of grand minima in 
the solar activity time-series, e.g. the Maunder minimum. 

These models are of course accompanied by amplitude varia- 
tions of the differential rotation. It will therefore be extremely in- 
teresting to measure the differential rotation at various instances 
of the 400-year record of telescopic sunspot observations. 

An attempt to determine the differential rotation during the 
period of the Maunder minimum indeed indicated a differential 
rotation that was stronger in the period of 1701-17 19 than on the 
present Sun (Ribes & Nesme-Ribes 1993). The whole rotation 
period of the Sun appeared to be shorter, with a slow-down at 
the turn of the 17th century. 

We investigated the period of 1749-1799 with the attempt 
to determine the differential rotation of the Sun in a thorough 
statistical analysis. 

2. Data set 

We used the drawings of the full solar disk made by Johann 
Caspar Staudacher in Nuremberg in the second half of the 18th 
century as described by Arlt (2008). All drawings were made 
employing the projection method which can be deduced from 
the orientations of the drawings of solar eclipses. The drawings 
lack any indication of the orientation. The sunspot positions de- 
termined by Arlt (2009) were obtained by matching the solar 
rotation profile with the spots in adjacent drawings. Here we 
cannot take the sunspot positions obtained by Arlt (2009) into 
account, since these employed the rotation profile by Balthasar 
et al. (1986) to find the position angles of the solar disks, and do 
not allow an independent determination of the differential rota- 
tion. 

The entire set of sunspot drawings by Staudacher was exam- 
ined in suitable pairs, which contain several apparently identical 
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Fig. 1. Association of spots in the two superimposed images of 
1761 June 23 (spots coloured in green) and 1761 June 24 (spots 
coloured in red). The difference in position angles determined 
by rotational matching is 29.8°. 



spots in both images. It is obvious that we have to compile a 
considerable number N of pairs in a sample to determine the 
differential rotation. Apart from the intrinsic plotting errors in 
the drawings, additional errors may occur when spots in adja- 
cent drawings were associated with each other, while they were 
actually not related. 



3. Data analysis 

3. 1. Defining the unknowns 

The entire 999 individual drawings of the solar disk were in- 
spected manually for suitable pairs of drawings. Drawings with 
less than three spots were omitted, since they do not deliver 
enough information to statistically determine of the free param- 
eters. A final set of 288 pairs were found with at least three spots 
in common. Several drawings were used twice if the combina- 
tions with the preceding and the succeeding drawings provided 
suitable pairs. As in Arlt (2009), the drawings were cut out of 
the photos by defining the left, right, lower, and upper edge of 
the solar disk on the screen. Parts of any elliptical deformations 
can be compensated for in this way. The actual spot pairs were 
picked in an image that contained the two superimposed draw- 
ings, which were now perfectly circular and slightly magnified. 
The scale in the disk centre is 0.12° per pixel in heliographic 
coordinates. Using the rotation profile by Balthasar et al. (1986), 
we estimated the expected displacements of spots. After one day 
a spot at 10° heliographic latitude becomes displaced by 0.63° 
in longitude against a spot at 30° latitude. This displacement cor- 
responds to about 5 pixels and should be measurable. Obviously, 
it is difficult to detect the differential rotation signal in the ob- 
servations because of these quite small displacements. It is the 
statistical multitude of drawings that delivers a measurement of 
the differential rotation with reasonable accuracy. 



Table 1. Comparison of the number of measurements versus the 
number of free parameters for given spot numbers and differ- 
ent methods using pairs of drawings. The maximum number of 
common spots that appeared in Staudacher's drawings is 20. 



Spots 


Measurements 


Parameters 


2 


8 


8 


3 


12 


10 


4 


16 


12 


5 


20 


14 


6 


24 


16 


7 


28 


18 


20 


80 


44 



In any i-th pair of drawings, we measured the cartesian co- 
ordinates (x\j, y\j) and (x\ 2 j , Vij) °f n i s P°ts visible in both 
drawing (1) and drawing (2), normalised in the solar radius, i.e. 
taking values between —1 and 1. We assumed that the longitudes 
Xij, the latitudes bij, and the two position angles of the draw- 
ings pi, qi are unknown, where j counts from one to the number 
of common spots rii. We also assumed that the profile of the 
angular velocity can be expressed in terms of 

o = n cq + sn sin 2 b. (i) 

In order to place a given spot at the same longitude in a pre- 
defined rotating frame, we need to know the precise times of 
the two drawings. Unfortunately, these are not always provided. 
Even when there are notes about the times, they are hardly more 
precise than to the full hour. In principle, we can also keep the 
time difference At, between the two drawings unknown. Note, 
however, that the three quantities f2 e q,i> Sili, and Ati are highly 
correlated. Any change in Ati can always be compensated for 
by a corresponding change in f2 e q,i- 

In the same way, if the spots are on similar latitudes, the 
angular velocity f2 e q,j and the differential rotation parameter 5Qi 
are highly correlated, since a lower equatorial velocity can be 
compensated for by a larger differential rotation parameter. We 
have a few options to reduce this redundancy: 

- We may assume that the times given by Staudacher are suf- 
ficiently accurate and fixe Ati. The parameters Q eq ,i, 5£li 
have to reproduce the longitudinal motion of the spots alone. 

- We may assume that the equatorial rotation period fi eq of the 
Sun has not changed and fix it at today's value. If this cannot 
explain the spots' longitudinal displacement, At, will com- 
pensate for this. Remaining discrepancies can be mediated 
by <5Qj. 

- We may assume that the total angular momentum of the Sun 
has not changed since Staudacher. If we assume that the 
shape of the internal rotation profile has not changed (except 
its amplitude), and if we assume that we know this profile 
precisely, we can find a relation for <5f2j(f2 e q,j)- Hence, only 
Q eq j and Ati are the unknowns. 

The first option is the most natural choice but turned out to 
be problematic, since the times given by Staudacher appear to 
be rough estimates only. There is also a considerable number of 
drawings for which no time is given except the date. Excluding 
these would leave us with too few pairs to obtain a meaning- 
ful result. The last option appears to be the most physical one 
since it is highly unlikely that the total angular momentum has 
changed over 250 yr. The solar wind carries away only about 



2 



R. Arlt and H.-E. Frohlich: Solar differential rotation in the 18th century 



10 30 gem 2 s~ 2 sr" 1 (Pizzo et al. 1983). In 250 yr this is less 
than 10~ 8 Aq, where A@ is the total angular momentum of the 
Sun of about A Q = 2 ■ 10 48 gem 2 s" 1 (e.g. Komm et al. 2003). 
The problem with this approach, however, is that fixing the to- 
tal angular momentum does not fix the distribution of fl at the 
surface well enough because of the low density and the less well 
known internal rotation profile of the Sun, which enters the total 
angular momentum. This is the reason why we chose a compro- 
mise of keeping the observing time open (Aij is a free param- 
eter) but fixed the rotation period at the equator to the present 
value. Because the specific angular momentum depends on the 
distance of a fluid element from the rotation axis, a change of 
ft(b) at a given total angular momentum means a considerable 
change at high latitudes, but only a small change at the equator. 
Assuming a constant 57 oq is therefore a fair approximation of a 
constant total angular momentum. 

From the two drawings of each pair and the two measured 
coordinates, we have Ani measurements for any given pair. The 
heliographic positions of these spots, the differential rotation 
parameter, and the time difference lead to 2rii + 4 unknowns. 
TableQ]gives an impression of the number of measurements ver- 
sus the number of free parameters. 



3 or more spots 
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Fig. 2. Marginal distribution for the differential rotation parame- 
ter for a minimum of three spots. A total of 156 pairs were used. 
The average value is 5Cl = — 0.048±g;g|| d _1 . 



4 or more spots 



3.2. Bayesian inference 

We used a Bayesian approach to determine the differential rota- 
tion. The method is based on the studies by Frohlich et al. (2009) 
and Frasca et al. (201 1), who determined the differential rotation 
of stars from photometric data. The idea is to define a model with 
free parameters and to determine the likelihood for many sets of 
parameters. The probability for each realisation of the model to 
have produced the data is computed. There is yet another free- 
dom in the system, which describes the uncertainty of the mea- 
surements. We have set this to a fixed value of 0.05 solar radii 
for the standard deviation, assuming Gaussian errors in all our 
runs. It is to represent the plotting accuracy when the spots are 
drawn using the projection method. 

A main concern of Bayesian parameter estimations is to ob- 
tain a comprehensive idea of the probability distribution in the 
whole parameter space. We therefore used a likelihood function 
for the i-th pair of drawings, 

A(a; i! y J ;A i , bj,Sn, At) = 



TT / — ex p 



( [xj - f x (x 3 ,b 3l sn,At)Y 

I 2a 2 

[yj- fv(^,b 3 ,sn,At)f 

2a 2 



(2) 



(omitting the indices i for the sake of legibility) where f x and f y 
compute the theoretical cartesian coordinates of all m common 
spots in the normalised solar disk from the set of parameters. In 
principle, we seek the probability distribution of a Stt common 
to all pairs of drawings in the investigated period. This problem 
would create a huge parameter space with the entire heliographic 
spot coordinates, all position angles and all time offsets plus a 
single parameter 5il. Three spots in three images would then 
deliver 18 measurements versus 12 unknowns. 

Such a problem is computationally extremely demanding. 
Fortunately, the problem splits into subsets of parameters. These 
subsets would be completely independent if the spots in one 
pair are different from the spots in the following pair. Then the 
probability distributions of the two individual parameters scans 




Fig. 3. Marginal distribution for the differential rotation parame- 
ter for a minimum of four spots. A total of 1 19 pairs were used. 
The average value is dil 



for the two pairs can simply be multiplied. In practice a num- 
ber of spots were "reused" in a following pair and thus lead to 
dependence on the previous pair. In principle, we could have 
used longer "chains" of observations. However, because sunspot 
groups evolve, it becomes arbitrary to associate common spots 
over more than two drawings. Problems may be introduced by 
erroneously marking spots as being identical in different draw- 
ings that are in fact not. We therefore decided to restrict the anal- 
ysis to pairs of drawings and to use the product of all probability 
distributions obtained from them for the final estimate of Stt. 

All possible heliographic longitudes from — tt to ir were 
taken into account as well as all possible latitudes, but in the 
sin b space to account for the smaller areas covered by given lat- 
itude bands at high latitudes. The range of 5il runs from —0.5 
d _1 to 0.5 d _1 and includes an increase of angular velocity with 
latitude as well as even a pole rotating in the opposite direction 
from the equator. While the latter situation is highly unlikely, 
we did not want to restrict the parameter space to certain solu- 
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3 or more spots 




-0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 

6(1 [rad/d] 



Fig. 4. Marginal distribution for the differential rotation param- 
eter for a minimum of three spots. A total of 55 pairs were used. 
The average value is SQ = -0.073±° °^ d _1 . 



tions. This is not a waste of computing power since the Markov 
chains will practically never enter this region of the parameter 
space because of the extremely small likelihood. The time offset 
is bound to —0.6 d < Atj < 0.6 d. The lower limit represents 
the situation of Staudacher observing near sunset on the first day 
and near sunrise on the second day. Formally, this corresponds 
to observing times at 19 h 12 m on the first day and at 4 h 48 m on 
the second day, assuming symmetry around midnight. The upper 
limit for Ati is the opposite extreme case with an observation at 
4 h 48 m on the first day and at 19 h 12 m on the second day. The 
extreme sunrises and sunsets at Nuremberg are at 3 h 52 m and 
20 h 10 m , respectively. These are not fully covered by our Ati 
range but it would only be relevant if both extremes were used 
for two consecutive observations and only during a few weeks 
during the summer. 

The full scan of a parameter space with more than a few 
dimensions is technically not possible. An efficient tool for ob- 
taining a representative idea of the probability distribution are 
Monte-Carlo Markov chains (MCMC, cf. Press et al. 2007). The 
enormous advantage of such a search is that one obtains confi- 
dence intervals for all parameters within the assumption of the 
model. Typically 64 Markov chains are "exploring" the parame- 
ter space in 20 million steps. In order to enhance the efficiency, 
MCMC was performed in an orthogonalized space, in which the 
new variables, a linear combination of the original parameters, 
were decorrelated. A principal component analysis (PCA) was 
performed using singular value decomposition (SVD, cf. Press 
et al. 2007) at several instances during the Markov chain walks. 
While the classical PCA relies on the covariance matrix and pro- 
vides as eigenvalues the square of the standard deviations, the 
SVD results in standard deviations themselves. An SVD is to be 
preferred because it is numerically more robust and faster (see 
e.g. Wall et al. 2003 for more details). The background math- 
ematics behind the SVD method used is very complex and be- 
yond the scope of this paper. Note that a full decorrelation can 
only be obtained if the probability distribution is a multivariate 
Gaussian; for all other distributions, SVD (or PCA) does reduce 
correlations between parameters, but does not make them vanish 
completely. 



3 or more spots 
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Fig. 5. Marginal distribution for the differential rotation parame- 
ter for a minimum of three spots. A total of 101 pairs were used. 
The average value is Sn = -0.039^°;°,^ d -1 . 



The set of parameters finally adopted was derived from the 
marginal distribution, which is the probability distribution of one 
parameter when integrated over all other parameters. Because 
we obtained the probability distributions from N pairs of draw- 
ings, we multiplied all distributions to obtain a total distribution 
for the differential rotation parameter of the Sun. The result tells 
us for a given model the probability of each configuration of this 
model to have produced the data. 

4. Results and discussion 

The probability distributions of all the 288 pairs were computed 
based on almost 1.3 • 10 9 Markov chain steps for each of the 
pairs. Several distributions show peculiarities, such as a maxi- 
mum at one of the edges of the considered parameter range. The 
reason for this is the association of spots in two drawings that 
have actually nothing in common. To exclude these unsuitable 
distributions, we selected only those for which the maximum 
probability lies inside the the tested intervals for Sil and At and 
not on the boundaries set in Sect. 13.21 We also excluded totally 
flat probability distributions by setting a maximum standard de- 
viation of a (supposed) Gaussian as a selection criterion. 

Figure |2] shows the marginal probability distribution for all 
relevant 5Q to have produced the data as derived from pairs with 
at least three common spots. The expectation value and the ap- 
proximate la confidence interval are Sft = — 0.048±0.025d _1 . 
If we restrict the data to pairs with at least four common spots, 
we obtain the distribution shown in Fig. [3] The expectation value 
is dn = -0.045 ± 0.026 d" 1 . The width of the distribution is 
slightly broader because fewer pairs define the result. However, 
since these combinations fix the unknowns better because of the 
larger number of spots, the increase in width is marginal. 

The selection of pairs with five or more common spots de- 
livers a very wide marginal distribution with an uncertainty of 
more than 0.05 in 5il. The number of pairs used is then 86. 

The results are compatible with the differential rotation of 
the present Sun of 5fl = —0.0501 d _1 for a sin 2 -profile, ac- 
cording to Balthasar et al. (1986). The deviations of the expecta- 
tion values based on Staudacher's observations from that value 
are not significant. Both probability distributions in Figs. [2] and [3] 
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Fig. 6. Butterfly diagram from 1512 spot latitudes determined through the process of Bayesian inference of the differential rotation. 
Individual spots were broadened as if they had persisted for 80 days to obtain a smoother picture. 



show a slight tendency to a smaller differential rotation in 1749- 
1799 than today. This may be worth keeping in mind if more 
historical data are recovered, which would improve the statistics 
of the present work. Even if our study does not find a signifi- 
cantly different solar rotation profile, it is remarkable that these 
historical sunspot observations made by an amateur astronomer 
with a rather simple telescope reproduce the solar differential 
rotation nicely. The result is an indication for the quality of the 
drawings, which also led to the conclusion of an unusal butterfly 
diagram by Arlt (2009). 

It may be worthwhile to check whether the average from the 
posterior distribution is a good measure of location for the differ- 
ential rotation parameter. We find from simulated distributions of 
known differential rotation that the average is indeed very close 
to the true value if the posterior distribution becomes very small 
near the boundaries of the considered parameter interval. This 
is not necessary (and is in fact not true) for the individual dis- 
tributions we obtained from individual pairs of drawings. The 
average of their average differential rotation parameters is not 
equal to the average obtained from the total posterior, and thus 
not close to the true value. All total posteriors shown in Figs. [2^5] 
are extremely small close to the interval boundaries of —0.5 and 
0.5. 

An interesting experiment is the splitting of the set of ob- 
servations into two smaller sets. Figure |4] shows the marginal 
distribution for the drawings of 1749-1761, while Fig. [5]shows 
the distribution for 1762-1799. The differential rotation param- 
eters are SCI = -0.073^;^ d" 1 and dn = -0.039^;^ d _1 , 
respectively, and indicate a decrease of the absolute value of the 
differential rotation during the second half of the 18th century. 
The difference is hardly significant though, because it is a la- 
effect. 

A large differential rotation bears the possibility of a hydro- 
dynamic shear instability, as studied first by Watson (1981), who 
found the solar shear to be nearly marginal. A 3D study with 
viscosity by Arlt et al. (2005) led to higher limits for the dif- 



ferential rotation, while higher powers of sin b may alter these 
results (Dziembowski & Kosovichev 1987), which were not in- 
vestigated here. If more data will become available and lead to a 
more accurate Stt, a test on its hydrodynamic stability would be 
required. 

A change of differential rotation could also be important be- 
cause the butterfly diagram has a different appearance during 
cycles and 1 (Arlt 2009). The equator was highly populated 
by spots and there is no clear migration direction of spot emer- 
gence. We can now use the spot latitudes determined as a side- 
product of the differential rotation determination and plot a but- 
terfly diagram. Fig.|6]shows the positions of the 1512 spots used 
in this analysis. Since we can only plot the positions of spots 
associated with each other in pairs of drawings, their number is 
much smaller than in Arlt (2009). The main features of that ear- 
lier paper are also visible in the present graph: relatively normal 
cycles 3 and 4, while cycles and 1 show a highly populated 
equator, which persists also into cycle 2. The appearance of the 
first two cycles of Staudacher's period somewhat resembles a 
topology that dynamo theorists call symmetric solutions of the 
magnetic field, in contrast to the anti-symmetric situation that we 
observe today. Symmetry here is meant with respect to the equa- 
tor. Of course we do not have magnetic field polarities at hand 
for historical observations, and consequently the true topology 
of the dynamo mode will hardly ever be accessible. 

Mean-field models typically show that symmetric and anti- 
symmetric solutions have very similar excitation thresholds (e.g. 
Chatterjee et al. 2004; Dikpati et al. 2004; Bonanno et al. 2005) 
and it is not straightforward to assume the one should dominate 
the other all the time. A small change in <557 might have favoured 
the symmetric solution in the 18th century, which exhibited a 
different butterfly diagram. Mean-field dynamo models with a 
back-reaction on the differential rotation tend to show excursions 
to symmetric geometry for reduced differential rotation though 
(e.g. Tobias 1997; Kiiker et al. 1999; Bushby 2006). Before we 
speculate too much, we conclude here that there is an indication 
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for a change of Stt during the 18th century, and we aim to stim- 
ulate future studies on the possible competition of the individual 
dynamo modes. 
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